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Abstract 

We investigate bright matter-wave solitons in the presence of a spatiahy varying nonhnearity. 
It is demonstrated that a translation mode is excited due to the spatial inhomogeneity and its 
frequency is derived analytically and also studied numerically. Both cases of purely one-dimensional 
and "cigar-shaped" condensates are studied by means of different mean-field models, and the 
oscillation frequencies of the pertinent solitons are found and compared with the results obtained 
by the linear stability analysis. Numerical results are shown to be in very good agreement with 
the corresponding analytical predictions. 



1 



I. INTRODUCTION AND MODEL 



In a mean-field theoretical framework the macroscopic wavefunction of a dilute gaseous 
Bose-Einstein condensate (BEC) is governed by a classical nonlinear evolution equation, 
namely the Gross-Pitaevskii (GP) equation The nonlinearity in the CP model is intro- 
duced by the interatomic interactions, which are taken into regard through an effective mean 
field. The coefficient (coupling constant) g of the nonlinear term in the GP equation is con- 
trolled by the s-wave scattering length a, whose sign and magnitude determine many of the 
fundamental properties of Bose-Einstein condensates (BECs) such as their shape, collective 
excitations or statistical fluctuations [11. Importantly, the scattering length can be tuned by 

w n it 

an external magnetic optical p|, or dc-electric field. The possibility of controlling the 
interatomic interactions in BECs has inspired many experimental and theoretical studies. 
The former include (among others) the generation of bright matter- wave solitons j^, ^ , while 
the latter predict that time-dependent scattering lengths can be employed to arrest collapse 



a 



in two LZH and three [8[ dimensional attractive BECs or to create robust solitons [9] , periodic 
waves [10| and shock waves 

On the other hand, recently, the possibility of varying the scattering length spatially 
has been proposed, utilizing, e.g., an inhomogeneous external magnetic field in the vicinity 



of a Feshbach resonance 



12| . These, so-called, "collisionally inhomogeneous" BECs have 



recently attracted much attention, as they are relevant to many interesting applications such 
as adiabatic compression of matter-waves , Bloch oscillations of matter-wave solitons 

12 1, atomic soliton emission [l^, enhancement of transmittivity of matter- waves through 



barriers 
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16|, dynamical trapping of matter- wave solitons 



15j, and so on. Moreover, 



rigorous mathematical results concer ning the existence and the stability of solutions of the 
relevant GP equations also appeared 17 1. 

In this work, we consider an elongated attractive BEC, with inhomogeneous interatomic 
interactions, which is confined in a highly anisotropic harmonic trap, such that the conden- 
sate is confined solely in the transverse (r = a/x^ + y"^) direction and free in the longitudinal 
[z] one. Then, its longitudinal mean-field wavefunction (j){z, t) satisfies the following nor 

malized one-dimensional (ID) GP equation {12!. 
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idt(f){z,t) 



'^dl-g{zM{z,t)\' 



(1) 



where space, time and density are respectively measured in units of the transverse harmonic 



oscillator length a_L = ^/hjmuj^ {m is the atomic mass and uj± is the transverse confining 
frequency), the inverse frequency cjj^, and the inverse length (2|ao|)^^, where oq is the (neg- 
ative) scattering length of the corresponding collisionally homogeneous system. In this work, 
we assume a collisionally inhomogeneous condensate, such that the nonlinear coefficient g 
in Eq. ([1]) has a spatial dependence of the form, 

giz) = 1 - e^^ (2) 

in the region z G being a small parameter. Such a setting may be 

realized in a lithium condensate, by applying an external magnetic field with a dc value 
corresponding to the minimum of the scattering length of "^Li [sj], and a linear gradient. 



which is controlled by the parameter e (see 15j for a detailed discussion). Note that the 



number of atoms Af in the condensate is given by Af= (a_L/2|ao|)-^; where N = \(j)\'^dz 
is the integral of motion (norm) for the normalized GP Eq. ([T]). 

As is well known, in the homogeneous limit of e = 0, the GP equation becomes the com- 
pletely integrable nonlinear Schrodinger (NLS) equation, possessing bright soliton solutions. 
In the inhomogeneous case e 7^ 0, soliton solutions can still be found and studied analyti- 



cally (for sufficiently small e) by means of the adiabatic perturbation theory for solitons 18 1 
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15|- 



(see below). Importantly, in this adiabatic regime, bright solitons feel an effective trapping 
potential induced by the inhomogeneity, in which both the center and the amplitude of the 
solitons perform oscillations when displaced from z = 

The purpose of this work is to demonstrate that these persistent oscillations are associated 
to the existence of a discrete eigenvalue in the corresponding linear eigenvalue problem. This 
eigenvalue, as well as the associated eigenmode, which is actually the translation mode of 
the solitary wave, manifests itself due to the inhomogeneity-induced perturbation (see also 
[2^ for a relevant discussion concerning the so-called internal modes of the solitary waves 
in nearly- integrable systems). We will treat the eigenvalue problem analytically using an 
approach based on the general theory of perturbed Hamitlonian dynamical systems of the 



nonlinear Schrodinger type as developed in 



21 



22l | (see also references therein). This way. 



we will derive the above mentioned discrete eigenfrequency associated with the translation 
mode. Moreover, we will demonstrate that this eigenfrequency is, in fact, the strength of 
the effective harmonic trap induced by the inhomogeneity, and, as such, it can directly 
be derived employing the adiabatic perturbation theory for solitons. Numerical results are 
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found to be in excellent agreement with the analytical predictions. 

Finally, we also study a modified ID GP model, which takes into account the effect 
of dimensionality. In fact, we explore the so-called nonpolynomial Schrodinger equation 
(NPSE) j^], which can effectively describe the longitudinal wavefunction of a truly three- 
dimensional (3D) "cigar- shaped" BEG. The NPSE model can be expressed in the following 
dimensionless form: 



idt 



_\q2 _ 3c/(z 



2^^ 2(l-^?(^)|0P)i/2j-' 
in which space, time and density are normalized as in the ID GP Eq. ([1]), the number 
of atoms is again N'= (a_L/2|ao|)A^, while g{z) is given by Eq. ([2]). We consider the linear 
eigenvalue problem for the NPSE model as well, and show that the deviation from the purely 
ID regime results in an eigenfrequency upshift. The value of the pertinent eigenfrequency 
will be compared to the oscillation frequency of a bright soliton in the full 3D GP model. It 
is shown that the agreement between the two is fairly good for sufficiently small values of 
the normalized number of atoms N , i.e., sufficiently below the collapse threshold. 

The paper is organized as follows: in section II we present the analytical and numerical 
results pertaining to the ID GP equation. Then, in section III we study the effect of 
dimensionality on the translation mode's frequency. Finally, in section IV we summarize 
our findings. 

II. ONE-DIMENSIONAL CONDENSATES 

Let us follow the approach, based on general theory of perturbed Hamiltonian eigenvalue 



problems, of 



22| to show the existense of the translation mode of the collisionally- 
inhomogeneous matter-wave soliton, and calculate its eigenfrequency. First we note that in 
the unperturbed case of e = 0, the GP Eq. ([T]) possesses an exact stable stationary bright 
soliton solution of the form, 

0bs(^; t) = r]sech[r]{z - zq)] exp{-ifit), (4) 

where rj is the soliton's amplitude and inverse width, zq is the soliton center, and fi = 
— (l/2)?7^ is the soliton's chemical potential. In the case e 7^ 0, the integrability is broken 
but Eq. ([1]) is still a Hamiltonian system, with Hamiltonian H = Hq + eHi, where Hq and 
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Hi are given by [for g{z) given by Eq. ([2])]: 

/OO i poo 1 

-{\d,^\'-\^\^)dz, H, = J -z'\<p\'dz. (5) 

The condition for the solution of Eq. (jl]) to be sustained under the considered perturbation 
is that the solution remains an extremum of the perturbed energy Hi [21], which in our 
case happens for = 0. The stability of the perturbed solitary wave is then determined 
by the location of the eigenvalues associated with the translation and phase invariance 
(the symmetries of the unperturbed problem which control the near-zero eigenvalues of the 
linearized equations). Having in mind that the four relevant eigenfrequencies were located 
at the origin c<j = of the spectral plane of eigenfrequencies u = u^ + iui (in the unperturbed 



system), one may follow [2l|, |22| and find the new location of the eigenfrequencies (in the 




D=| ^"'^ I V |. (8) 



perturbed system) by means of the equation: 

det(eM-c<;2D) = 0, (6) 
where the matrices M and D are given by 



and 



-^0,3,2^^ 

In the above equations, star denotes complex conjugate, <, > is the inner product, and 
6Hi/6(j)*' is the functional (Frechet) derivative. Note that the matrix M is generally diagonal, 
with nonzero elements mu and 17122 acounting, respectively, for the perturbation-induced 
breaking of the translational and phase invariance. However, in our case, the considered 
form of the inhomogeneous nonlinearity does not break the phase invariance of the system 
and, as a result, 17122 = 0. As we will show below, a consequence of the phase invariance of 
the system is the appearance of a double zero eigenfrequency in the linear spectrum. 

The nonzero elements of the matrices M, D can be directly calculated and the results 
are inn = (2/3)?7'^, du = i] and ^22 = —i]^^- Thus, Eq. ([6]) leads to the following simple 
algebraic equation. 



Equation (Q provides the double zero eigenfrequency cu^ = reflecting the phase invariance 
of the system, as well as the new location of the eigenfrequencies which were associated to 
the translational invariance (which is now broken due to the presence of the spatially inho- 
mogeneous nonlinearity) . In fact, the latter pair of eigenfrequencies provides the frequency 
of the translation mode (which will be called Q) which is given by Q"^ = {2/?))trf. To this 
end, using the relation /i = —{\/2)rf, the latter result can be expressed as 



^ = y-^e/i. (10) 

It is interesting to note that an alternative analytical approach can be used to obtain 
the value of the eigenfrequency. This approach is based on the adiabatic perturbation 
theory for solitons ^Jj], which states that approximate soliton solutions, characterized by 



parameters that are unknown functions of time, can stil 
([1]). Following the methodology expounded in Refs. jl2. 



be found for the perturbed system 
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15| , it is straightforward to find 



that the soliton's center evolves according to the following equation of motion, 

where 77(0) is the initial soliton amplitude. Taking into regard that = — (l/2)?7^(0), it is 
readily found that Eq. flTTl) becomes (Pz^/dt^ = —Q'^zq, where Q is given by Eq. flTU]) . This 
means that the soliton center behaves like a Newtonian unit-mass particle in the presence 
of the effective trapping potential Vcs = {1/2)Q'^Zq. Thus, when displaced from the trap's 
center (zq = 0), the soliton will perform harmonic oscillations with frequency Q, which is 
nothing but the eigenfrequency of the translation mode of the solitary wave. This result 
indicates the physical significance of the translation mode's eigenfrequency, which is the 
same as the soliton oscillation frequency in the effective trapping potential induced by the 
inhomogeneous interactions. 

The above analytical predictions have been checked by two different types of numerical 
simulations, in which Q was directly derived by a linear stability analysis of Eq. ([1]), or 
obtained as the oscillation frequency of the bright soliton in the framework of the GP model 
ofEq. (HD. 

Let us first discuss the results obtained by the linear stability analysis, which can be 
performed upon considering small perturbations around the unperturbed soliton of the form 

0(z, t) = [0bs(^) + e (M(z)e-*-* + i;*(^)e*"**)] e"*'^*, (12) 
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where u and v represent the normal modes oscillating at eigenfrequencies ±00. Substituting 
Eq. f|T2l) into Eq. ([1]), we obtain the following Bogoliubov-de Gennes (BdG) equations (valid 
to leading order in the small parameter e): 



where = 2g{z)(f)lg and F2(0^s) = fl'(^)0bs' ^^^^ solitary wave solutions 0bs- 

The above BdG equations have been solved numerically to find the eigenfrequencies, and the 
resulting spectral plane Re(co')-Im(c(j) is shown in the left panel of Fig. [T]for e = 0.01, = 
1.2, and fi = —0.18. Note that the eigenfrequencies appear in pairs due to the Hamiltonian 
nature of the system under consideration. As shown in the left panel of Fig. [U there exists 
a pair of eigenfrequencies at the origin (corresponding to the symmetry associated with 
the phase invariance), as well as a pair of eigenfrequencies located at ±Q, i.e., in the gap 
between zero and the continuous branch > — /i. The latter pair is the translation mode's 
eigenfrequency appearing due to the presence of the inhomogeneous interactions that break 
the translational invariance of the system. 

Moreover, in the right panel of Fig. [1] (see solid line) we plot f2 as a function of e, for 
= 2.8 and fi = —1. It is clear that the eigenfrequency Q obtained by solving the BdG 
equations perfectly follows the square-root law of Eq. (fTOj) (as well as the prediction of the 
perturbation theory for solitons); in fact, the respective curves are identical and cannot be 
distinguished from each other. Note that in this figure e G (0,0.16) so that the spatially 
inhomogeneous nonlinearity remains attractive (see also the discussion below). 

Next, we have numerically integrated the GP Eq. ([T]) with the initial condition found 
using a fixed point algorithm (Newton- Raphson), with the initial guess being the soliton 
solution in Eq. (jlj) at t = 0. This "exact" stationary soliton solution was subsequently 
displaced so that zq 7^ 0. Note that for a fixed value of the chemical potential fi (or soliton 
width r]~^), the displacement of the soliton is such that the soliton oscillates in the region 
where the nonlinearity is attractive, i.e., in the interval z G {— e~^/^, e"^/^}. An example of 
the initialization of the system is shown in Fig. [21 where both the inhomogeneous nonlinear 
coefficient g{z) and the density of a bright matter- wave soliton (displaced from the effective 
trap center) are shown. The resulting frequencies of the soliton oscillations are depicted 
by stars in the right panel of Fig. [H which are clearly located very close to the solid 



—uv = 



UJU = 




(13) 



(14) 
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FIG. 1: (Color online) Left panel: The linear spectrum of a bright soliton of the GP Eq. ([T|) for 
e = 0.01, N = 1.2, and fi = —0.18. The bright soliton's translation mode has an eigenfrequency 
= 0.049 in excellent agreement with the prediction of Eq. (jlOp . Right panel: Solid line shows 
the translation mode's frequency as a function of e obtained by the BdG Eqs. (jl3p - (jl4p . as well as 
the theoretical prediction of Eq. (jlOp for = 2.8 and fi = —1; the two results are identical and 
the respective curves cannot be distinguished from each other. Stars depict the soliton oscillation 
frequency obtained by the numerical integration of Eq. ([T]) . 

line Q{e). We conclude that there is a remarkable agreement between the solution of the 
BdG equations, the integration of the GP equation and the predictions of the two different 
perturbative approaches. 

As far as the validity of our predictions is concerned, we note the following: if the soliton 
width, r]~^, is sufficiently smaller than the characteristic width of the inhomogeneity, e~^/^, 
the soliton satisfies the relevant predictions very accurately in the small and intermediate 
oscillation amplitude regime. This is expected to occur due to the robustness of the soliton, 
captured by the perturbation theory for solitons, in the perturbed-inhomogeneous-system, 
and the validity of the (linear) BdG analysis for small- amplitude oscillations. However, in 
the case t]^^ ~ e^^^^, or/and for large amplitude oscillations, the soliton evolves under a 
strong inhomogeneous perturbation and, as a result, nonlinear effects, as well as emission 
of radiation, become important (see, e.g.. Fig. 2 in Jisl). In such cases, our assumptions 
cease to be valid and, as a result, our analytical (perturbative) approaches should not be 
expected to agree with the numerical results. 
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FIG. 2: (Color online) Solid line shows the density of a bright matter-wave soliton with = 1.2 
and n = —0.18 (the parameter values are the ones used in Fig. [1]); the soliton is initially placed at 
zq = 2. Dashed line shows the nonlinear coefficient g{z) for e = 0.01. The initial displacement of 
the soliton, zq, which sets the amplitude of the soliton oscillation, is such that the inhomogeneous 
nonlinearity remains attractive in the region where the soliton oscillation takes place. 

III. CIGAR-SHAPED CONDENSATES 

In many experimentally relevant situations, the transverse confinement of the conden- 
sates is not sufficiently tight and, as a result, deviations from ID are quite relevant. In 
such cases, the condensates are "cigar-shaped" and their mean-field descri ptio n requires the 



consideration of either the 3D GP equation, or other effectively ID models 23|, |2J, |25| . Here, 



we consider the NPSE model of Eq. ([3]), which has successfully been used to describe recent 



experimental results 26 1. 



Following the same procedure as in the case of the GP Eq. ([T]), we introduce the ansatz 
of Eq. f|T2|) in Eq. ([3]) to obtain BdG equations similar to Eqs. f|T3l) - f|T4l) . but with the 
functions Fi and F2 given by: 

p ..2 ^_ 9g\z)4>l-Ug{z)<Pl + 4 

^i(<^bJ - 4[1-,(.)0L]^/^ ' ^''^ 

^^(^J = 4[i-,(.)0L]^/^ • ^''^ 

The eigenfrequency spectrum of the NPSE model has a form similar to the one pertaining 
to the GP model. However, as shown in the left panel of Fig. [3l the translation mode's 
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FIG. 3: (Color online) Left panel: Similar to the left panel of Fig. [H but for the NPSE model of 
Eq. ([3]). The inhomogeneity parameter is the same (e = 0.01) and the chemical potential is chosen 
(/U = 0.762) so that the norm of the bright soliton be the same in both GP and NPSE models. 
The translation mode's frequency is up-shifted, i.e., \Q\ = 0.062, as compared to the respective 
value calculated with ID GP equation (depicted by stars). Right panel: The translation mode's 
frequency as a function of the norm A^, calculated by the BdG analysis of the NPSE model 
(dashed line) and the ID GP model (solid line). The eigenfrequency O is clearly upshifted due to 
increase of the dimensionality. The inhomogeneity parameter is again e = 0.01. 

eigenfrequency is clearly upshifted: For the same value of the inhomogeneity parameter 
(e = 0.01), and the same norm (A^ = 1.2), we find that Q = 0.062, which is ~ 27% upshifted 
as compared to the corresponding value obtained for the GPE model {Q = 0.049). 

We have found that this effect, i.e., the upshift of the translation mode's eigenfrequency 
due to the increase of dimensionality, is a generic feature of the system, as shown in the 
right panel of Fig. [3l In this figure, the eigenfrequency Q is given as a function of the 
norm (for e = 0.01) for both the NPSE model (dashed line) and the ID GP model (solid 
line); the latter is apparently a near-straight line [due to the validity of Eq. ffTTl) and the 
fact that N = 21] for the bright soliton]. The range of values of N is such that N < 1.33, 
i.e., below the collapse threshold Nc = 1.33 {2^, and N > 0.7, so that the inhomogeneous 
nonlinearity remains attractive in the region where the soliton oscillation takes place. The 
latter requirement is due to the fact that is inversely proportional to the soliton width 
(see also Fig. [2]). 

To corroborate the validity of the NPSE analysis for the higher dimensional case, we have 
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numerically integrated the pertinent 3D GP equation 



--V^ + -r^ - q(z 
2 2 



^, (17) 



in which = r^^dr{rdr) + d'^, and the additional trapping potential term (l/2)r^?/^ in 
the transverse direction has been incorporated. The initial condition was obtained by a 
relaxation method, using the initial guess 



V2 



a 



r 



2 



^{t = 0) = ^exp ( -— ) 0bs(^) (18) 



where the width of the Gaussian in the transverse direction is a = a/1 + 1 0bs | ^ (as per the 
231]) and </)bs is the soliton profile of Eq. (jlj). An example of the simulations 



ansatz used in 



is shown in Fig. HI where the spatio-temporal plot of the soliton's longitudinal density is 
shown (for = 0.757 and /i = 0.92) and compared to the prediction of the BdG analysis 
of the NPSE model (dashed line). It can be observed that the agreement between the two 
is very good: The BdG analysis predicts an translation mode of frequency Vt = 0.033, while 
the result of the 3D simulation shows that the oscillation frequency is 0.03307, with the 
error being 0.2%). We note in passing that in the purely ID regime, the BdG analysis of 
the ID GP model predicts an eigenfrequency Q = 0.02989, or 10% discrepancy from the 
3D result in this case. We should also remark that the deviation of the ID NPSE result 
for the oscillation frequency from the fully 3D GP equation becomes worse for larger values 
of the norm A^, i.e., when approaching the collapse threshold of A^ = 1.33. For example, a 
discrepancy of 8% was found for A^ = 1 (the BdG analysis of the NPSE predicted Q = 0.0457 
while the 3D GP model provided an oscillation frequency 0.0498). 



IV. BRIEF SUMMARY 



We have investigated the static and dynamic properties of matter-wave bright solitons in 
a parabolic coUisionally inhomogeneous environment. It is shown that a translation mode 
of the soliton can be excited in this setting and we have found analytically its frequency. 
Numerical results, based on the relevant Bogoliubov-de Gennes equations, were found to be 
in excellent agreement with the analytical predictions. Moreover, we have shown that in 
the purely ID setting, the oscillation frequency obtained in the framework of the mean-field 
Gross-Pitaevskii model coincides with the translation mode's frequency. Deviations from ID 
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FIG. 4: (Color online) Spatiotemporal contour plot of the soliton's longitudinal density for N = 
0.757, fi = 0.92, and initial soliton position zq = 2. The dashed line shows the harmonic oscillation 
of the soliton center, with a frequency predicted by the BdG analysis of the NPSE model. Note 
that for the chosen value of e the inhomogeneous nonlinearity remains attractive for — 10 < x < 10. 

were also considered upon analyzing the nonpolynomial Schrodinger model and it was found 
that the translation mode's frequency is upshifted. Numerical results obtained by the 3D 
Gross-Pitaevskii equation demonstrated that the oscillation frequency of the bright soliton 
in this "cigar-shaped" setting is generally underestimated by the nonpolynomial Schrodinger 
model, which, in turn, is underestimated by the ID Gross-Pitaevskii equation for the same 
total number of particles. However, for numbers of atoms sufficiently below the collapse 
threshold the BdG analysis of the nonpolynomial Schrodinger model accurately predicts the 
soliton's oscillation frequency in the 3D setting. 

This work has been partially supported from "A.S. Onasis" Public Benefit Foundation 
(G.T.) and the Special Research Account of the University of Athens (P.N., G.T. and D.J.F.). 
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